{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 100000\n",
    "\n",
    "x = np.random.normal(size=m)\n",
    "X = x.reshape(-1,1)\n",
    "y = 4.*x + 3. + np.random.normal(0,3,size=m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXkAAAD4CAYAAAAJmJb0AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjMsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+AADFEAAAgAElEQVR4nO3dcYwc91UH8O+79TjZMyVrK0eJN77YlMqmrmtbOVKDkajTgkOTJm7S1C0JKgLJQiqiNsb03EZKgkA56URTpBYhCxBICek5OGycuuCmtREQ4dC73rnGjQ0pdS5eB2rqbFrsrb239/hjd9azuzOzM7szO7Mz348Uxbdzd/vz+vbdb97v/d5PVBVERJRMQ1EPgIiIwsMgT0SUYAzyREQJxiBPRJRgDPJERAm2JOoBWN188826evXqqIdBRDRQZmZm/ldVR+yuxSrIr169GtPT01EPg4hooIjIa07XmK4hIkowBnkiogRjkCciSjAGeSKiBGOQJyJKsFhV1xBR+hRmi5g8ehYXSmWszGWxb/ta7Nicj3pYicEgT0SRKcwWsf+5UyhXqgCAYqmM/c+dAgAG+oAwyBNRZCaPnm0EeFO5UsXk0bN9CfJpuItgkCeiyFwolX09HqS03EVw4ZWIIrMyl/X1eJDc7iKShEGeiCKzb/taZI1M02NZI4N929eG/txR3kX0E9M1RBQZMy0SVl7cLee+MpdF0Sag9+Muop8Y5IkoEN0uYu7YnA8lB94p575v+9qm64C3u4hBW6xlkCeinkW5iOkUdDtV7nRzFzGIi7UM8kTUs6hKId2CrlNu3Zqi8XsXEXXJZze48EpEPYtqEdMp6O49eBK5YcP2awS1Xw5+FWaLtjl8IN6LtQzyRNSzqEohnYJrVRX/96MF22sK+C6TNO8YnMR5sZbpGiLqmddFzCAXLQuzRQyJoKpqe72yaP844H3mbY7XaQYP9K/ks1sM8kTUM7dFTKdA2bpoaf0FkDWGcKWy2PjcZUsz+KMPb2g8jzmzdgrwnbjNvK3jFdRm/m6euH9DbPPxAIM8EQXEbhGzdWG0lXWHqfXzrAEeAC5fq2Lvsycbz2OXi/fKbebdOt5OAT4j0tUY+olBnogC9UjhFJ55+XXPs+xiqYy9B092/PzqojaqWHpZ6HSbefv95VFVZQklEaXHI4VTeOrEvK+vEcDzL4QLpXLHXLybXLZWcbPp8a+iVK4AAJYPG3j0Q+u7/uUR9xJKBnkiatLL4ugzL7/u+/n8hOqsMdRTLv6HVxfwu1NzsCaD3rxSwb6/raWCnFoddMISSiIaCGZOulgqQ3F9cdRrXXm3wderK5XFrnPxQC3ls2jzeKVaSwXZNUwzs+75XLZxJ9AqziWUPQd5EVklIsdF5BUROS0in6o/vkJEXhSR/6z/f3nvwyWiMPXafncQFiKdFEtlTB49iwduzyOfy0JQC+xP7tyEcxN346XxO/HYvesj65rZrSDSNQsA9qrqN0XkbQBmRORFAL8O4OuqOiEi4wDGAXw6gOcjogBZ0zNO83Cv6YiPv3eV75x8GARAbtjAm1cqvr6uWCrj0EzRcXE27K6ZYRAN+PZKRJ4H8IX6f+9T1TdE5BYA/6iqrr/uxsbGdHp6OtDxEJGzTiWOpnwui5fG7+z4vVrr4UWAJQJU7HIkIbOrcc8MCdQhZdMqlzUgApSuVGIfzEVkRlXH7K4FuvAqIqsBbAbwMoC3q+obAFAP9D/h8DW7AOwCgNHR0SCHQ0QdeCkZdEpHWO8AbsoauHxtAZVqc1hVBSrhpukd2T3tEICdW0bx5ZNvNKprnFivD0K3SSeBLbyKyI8BOARgt6r+wOvXqeoBVR1T1bGRkZGghkNEHnRKw+SyBm40hrBnag5bJ441FmBbF2hL5UpbgI+jyqLi+JmLmHv0l3Fu4m7kfSyYDurRgIGka0TEAPBlAEdV9XP1x86C6RqiWHAqi9w6ccyxZNDICKDNPWDMFEimyzr1uBDUKmK2rRvBoZmir4qdfC4bu3x8qOkaEREAfwHgFTPA1x0G8AkAE/X/P9/rcxHRdV7r2Z16rk+/dgmXr9p3agRgOzM3HxnkAA+gUR761Il5ZI0hLPexSGv+UhyUFE7PM3kR+QUA/wzgFNBYz/gMann5gwBGAcwDeFBVL7l9L87kibyxWzDNGhnbqhCn2bqX5ltpkTUyeOD2vO9ZPeBtUdpJUF05Q53Jq+q/4Pp+gVbv7/X7E6VBpzd76/XLVxc8n1DklHdngL+uXKni+JmLeOL+Ddg9Nefra7vd7dqvowS545UoYp12mdpdd6oMsZuxx3k3ZpwUS2Xs2Jx33NXqpNvXt9eNZ14xyBNFrNOb3W9nxM1/8NWmNgR2W/WpXUYEhdkiLl9zXqew0+1u134dmcggTxQxtwOnC7NF32/6N69U2vrN3Gjwrd5JVRW7p+b6VgraryMT+S9PFDG3N/X+5045HkjtxrwTMFM9frf3k3d7D57EmvEjTfsIvLC7wwqjDw6DPFHE3NIp5UoVqs6VDW4ulMp4/IXTPXVtpM6qql117NyxOY8n7t/Q1AwtjKMEA+9d0wuWUFLaeDkouttSR2Momp4xaddLSWW3+ta7hoi889oczEzn+D3MggE+GnE7QITpGqKI+GkOxgqZ+HHqnR+3klUGeaKIuM34WnO01vwtRS+fy+KPP7pxIA4QYbqGKCJu54nelDVw5doC9kzNNY6lM4M90L4DtptzSal7F+obp4D4HyDChVeiiHjNyZuyxhDKlUUMCWA2hswaQ6guKq4NQJvfuDEygp0/uwrHz1z0/UsyisVVN1x4JYoh60zQS5Ap11dSLZ1/G4+Rf5WqNo4qXLY0g8vXvP6ybU7JBNVkLCwM8kQhaz1ByTxSLjdsQBV4q8MJRRQ+rwE+3xLE+9VkrBcM8kQhag0C1sZi3IU6WHJZoy1F49Z3iEGeKCHcbte54zQ57O64+tVkrBcM8kRdKswW8djh044HPgOcrSeJXf27U2VTnGrlWSdP1AUzDWPX1928XR/EQ5/J2bZ1I22P9avJWC8Y5Im60Gm3arFUZu16whyaKbY1H+tXk7FeMF1D1IVOOVeen5o8Tguq1k1qccSZPFEX3HKubgGeb7j4MjLSsaVznBZUveLPHJEPhdkitk4cQ7FUtg0IQ9JhBt9NY3jqi0pVO56gFacFVa+YriHyqLXm3RrMzdi92CFH0+k6RcttB3HcFlS9YpAnauFU9+622MrYnWytO10HCYM8kYXbNvVBzMdS7wSIVTMyv5iTJ7Jw2qa+9+BJztYTbrnDgelDIr4O6I4bBnkiC6fZejVGLbkpePlcFo9+aL3t6VtVVV8HdMcN0zWUaq3592EfLWcpGcwFVTPfvvfgybZf6nFrOuYHZ/KUWmb+vVgqQ1HLvzPAp491h+qOzXksOty1DeqaDGfylAp2FTNeDtKmZFu2NNM2Ox+EpmN+cCZPiWc3Y98zNcfeMoTL16pYPX4EWyeONXLu29aNtO1ZG9QaeYAzeUoBuxk7l1HJyiyVnX7tEg7NFNs2uj1we7z707hhkKdEK8wWOWMnT8qVKp55+fW2RVcFcPzMxWgGFQCmayixzDQNkVdOpbKDuugKMMhTgnFhlfzKiH0HuUFddAUY5CmBrJ0iibzKGhl8/L2rYn/Sk1/MyVOiPFI4hadPzHNhlXzJZQ08du967Nicx9htKxwPZh9EgQR5EflLAPcA+J6qvrv+2AoAUwBWAzgH4KOq+mYQz0dkpzBbZICnriy7YUnThqhBDuqtgkrX/BWAu1oeGwfwdVV9J4Cv1z8mCs1jh08zwFNXBnlhtZNAZvKq+k8isrrl4fsAvK/+578G8I8APh3E8xEVZot47PBplMoVALWdi2xJQN0a5IXVTsLMyb9dVd8AAFV9Q0R+wu6TRGQXgF0AMDo6GuJwaNCZrQnsFlQZ4MnNEICbhg28eaXSdgbvoC+sdhL5wquqHgBwAADGxsZ4t022Wg/zIPLKuqgKOJ/8lVRhBvn/EZFb6rP4WwB8L8TnooRjzTv5JQCe3LmpLYAnbWG1kzCD/GEAnwAwUf//8yE+FyVA6wxr27oRHD9zERfqjcWI/FAgVcHcSVAllM+gtsh6s4icB/AoasH9oIj8JoB5AA8G8VyUTHZnqz51Yj7iUdEgE9R+rjoF+qSnb0RjdKzZ2NiYTk9PRz0MigB3qFIYMiJYVHUM3nZrPVkj03SQyCAQkRlVHbO7xrYGFAtJrlOm6FRVG2cI2J3T6nRw++TRs30cZbgY5ClyhdkihhwaQxEFxS54O00ukjTpYJCnSJm3y04tXomC1Bq8nTZBJWlzFIM8RYqlkRSErJHBw1tGkc9lIfDeMnjf9rWJ6zrZKvLNUJRObrtXibzI57KOFTFOC6qtwdv8miRX1zDIU6C8lKNx9yr1Kp/L4qXxOx2v+wneSd8cxSBPgbGrdTeP37O+iR5/4TQDPHXNazol6cHbK+bkKTBeytEKs0W8eaXS76FRggxaDXvUOJOnnljTM071MdaKhiTVH1P/+Sm0TfpOVq84k6eumemZYofeMorajtbCbDFR9cfUfwpvE4XWn02nzVBpwCBPXSnMFrH34EnPuXXzTXajwR856o2XiUIadrJ6xXQN+dbtBiYutlIQvGxUSsNOVq84rSLfuIGJouK1siYNO1m9YpAn39I4G6LwOS2qZkQgqNXGe62sScNOVq+YriHfVuay3KlKgcsNG/hRZTGQtr9p2MnqFfvJk2+F2SL2TM3xtCYKlAB4aMsonnn5dVRVkRHBlp9ajnPfL6c+UHfCfvIUqB2b83hoy2jUw6CEyQ0bODRTbCzoV1Xx0ncusQyyR0zXkGetm0uIuvH2ty3F9354relOUABPO6HNMkjO5r3jTJ48sdtcQtSNH/yoiofqbYGBWoD3k/rjwr8/DPLkCcsmKSjlShXHz1zES+N3Ip/L+l7b4V2kPwzy5AlnTxQk8+fJ789VWssge8EgT55w9kRBMn+e3H6u8rls02lPfurk6TouvFJHhdkiSleuRT0MSgjrbHzf9rW2JzgxmAeHQZ5c8RQnCoK5uJpvqXXnpqXwMciTI7PTpN9GZERWxhAw+eAmx8DNE5zCxSBPth4pnMLTJ+a5q5V6tiST6fxJFBoGeQLQvNEpN2zwiD4KDDcwRYtBPoVad65uWzeCQzPFRt6dAZ6CxhLc6DDIp0zrQmqxVGZahkLHEtzosE4+Zex2rjLAU1CMIYGRae4Mzw1M0eJMPoHcTqnnbTOFafLBjbX/syQyNhjkE8YuHbP/uVON6yIAKyIpDBmpzeBZEhkvDPIJ43RK/e8enMMigzuFqKramFAwyMcHc/IJUpgtOrYAZoCnfjDLJSk+GOQTwkzTEEWN6z7xEnq6RkTuAvAnADIA/lxVJ8J+zqRq3bCkCrxVrmBlLovLVxfYX4ZigeWS8RJqkBeRDIAvAvglAOcBfENEDqvqt8N83iRqXVC1bljiKU0UloyIr95FLJeMn7DTNXcAeFVV/0tVrwH4EoD7Qn7OROLJTNRv5ybuxh9/dCOyhnvvGbOqhv3e4ynsdE0ewOuWj88DeK/1E0RkF4BdADA6OhrycAYX85zUb6vHjyCfy+KB2/N46sS87ecIgO888cH+Dox8CXsmLzaPNd37qeoBVR1T1bGRkZGQhzO4vOQ5lw8bjcORrbOrXNYIdWyUXMVSGYdmilg+bP8zxPx7/IU9kz8PYJXl41sBXAj5ORPJ7gSdVqqw3V1YmC1i99Rc2EOkhCpXqrhhyRCyRqbtBCfm3+Mv7Jn8NwC8U0TWiMhSAB8DcDjk50ykHZvzeOL+DY4zKgAolSvY/9wpFGaLjcceKZzCHgZ46tFb5QqeuH8Dz1sdQKHO5FV1QUR+G8BR1Eoo/1JVT4f5nEn3o8qi6/VypYq9B09y5k6emSk+tyqtlbks2xUMqNDr5FX1KwC+EvbzJFFrozGvtfA8ro/8MFMuTulApmUGG3vXxJRdozGioC0fNppm55NHz6JYKjfq41sP3qbBwyAfU4+/cJp18RSqrJHBox9a3/iY6ZhkYu+aGCrMFn0dwTdUL1S1q1clcsKF03TgTD6GHn/BeW16+bCB4aVLbA9k+OnPfAULbDdJHuTrC6mUfAzyMfNI4ZTrLN68vTYXZM22rs9OzzPAkydcSE0XBvkYKcwWHbePA2jsXN337ElU6gG9WCrzQBDqKJ/L8ji+lGKQj5FOhy3cs/EW7Jmaazt4mwGeOnlp/M6oh0AR4cJrjHRqQnZoptgW4Ik64YJ8unEmH6HWzU43ZQ2Uys75eJZUUjc4MUg3BvmI2G12MjKCIWH6hYKVZ6fIVGO6JgKF2SL2HjzZNjOvVJWzLnLllnrJZY22Az5YSUMM8n1mzuCd+suw7Qw5yeeyeHLnJnx+5ybbYH7Pxltwo3H9LZ3LGtzwREzX9BuP8aNu5LJGW4WMdT1n27oRHJopNv1sXV1w71hK6cAg32c8xo+68VbLgnxrn5mtE8faJg/lShWTR89yJp9yDPJ9tjKXZUdJ8u2mrIGtE8ccNzQ5TR44qSDm5ANUmC1i68QxrBk/gq0Tx5pOaDLt276Wdcvk2+VrCyiWylDUKrFaTwBzOmuVZ7ASg3xAzAVVtzciULvNfmjLKAM9+VKpNq/Im6kY077ta1lZQ7YY5ANit6BqHsXXOrP/wx0b8OTOTVEMk2Ion8vi4S2jvr/OmooxzwDmGazUijn5gDjlPs1SyWKpjD1Tc5h+7RL+cEftzWeewkPpZc62zWD89In5pr0SWSODG40h286krakYHvpBdjiTD4iX3KcCeOrEfGNGz/x8urXWsY/dtgK5YaPt+qMfWs9UDHWNQT4gdjlRJ+ahIDs257nDNcWW3XD9Rtpc07HO2M06d6ZiqBdM1wTEfMOZG1SG6gch27G+kfMsqUwtc3EecF7TMevcmYqhbonGaB/92NiYTk9PRz2MQBRmi9g9Ned4PV+vdQbQ1KgMAJuUpYx5oIfdP7kA+O7E3f0eEg0YEZlR1TG7a0zXhGTH5jyyhvPLWyyVse/ZkwDQdCu+fNhggE+ZYn2Dkx3WuVOvGOS75GXj0xP3vwfGkPPSamVRsXtqDpNHz2Lf9rX47sTdbFCWQhkR1rlTaJiT74JdL3gzt2rNm1rz9G55d3NW/+z0vOuhITQ4li3NYOkS+9LHVlXVtjUdnsVKQWGQ70KnRTIr82O3/DxQm9W/9J1LwQ6U+i5rZBqVL62TASfmoR5cXKUwMF3TBb/NoMySSUoeY0iwfNiwLW1sLX1cPmy0pe+YkqGwcSbfBadOknaLZIXZoqdbdho8eQ8pldbZeeu5vkzJUNgY5Luwbd2I7fZzuxkZZ/HJ9Pmdm7oKzkzJUL8xXeNTYbaIQzPFpgAvAB643f7Ny1n84Nr6jhUwMu3VUQ9vGWWgpoHBmbxPdouuCuDLJ9/A8TMXeRueIOe+X8bkRzYyvUIDjTtefVozfsRTvxkjI1i2dAlLIgcYd5vSoHDb8cqZPPwthnk9vq9SVQb4AcfdppQEqc/Jez3RyeSn2yQNLpY2UlKkPsi7bWyyY9f2dbmlBzgNPrbypSTpKV0jIg8CeAzAzwC4Q1WnLdf2A/hNAFUAv6OqR3t5rrB42dhkl855afzOputedjZSNKz17I8UTuGpE/OOn3uOOXhKmF5n8v8O4H4A/2R9UETeBeBjANYDuAvAn4pILHMcnbr/eU3n3LAk9TdFsXWhVMbk0bMozBYxdtsKx9O48szBUwL1NJNX1VcAQKTtbXMfgC+p6lUA3xWRVwHcAeBfe3m+MOzbvrZtFp41Mti2bgRbJ47ZLrJa+9RwFh9/5i/nfX97ElDYVkcxB09JFVZ1TR7ACcvH5+uPtRGRXQB2AcDoqP8T63tl1/1v27oRHJopugZuM51jl9OneKpU7YtfMyLMwVNidQzyIvI1AD9pc+mzqvq805fZPGb7DlPVAwAOALU6+U7jCUPrVvOtE8c6Bm4zneOU06fBsWhp9UuUNB2DvKp+oIvvex7AKsvHtwK40MX3iYSXwG3e2nutm6f+8nOEIuvhKcnCWi08DOBjInKDiKwB8E4A/xbScwWu05s+lzUaMz/WzcdTRqSt74yREbb6pdTpKciLyIdF5DyAnwNwRESOAoCqngZwEMC3AfwDgE+q6sAkrvdtX+tYgSEAHrt3fePj1rr5XNbAsqXuQZ919eGrLCqWLV3StJ9h8iMbMfngxqbHmIunpGPvGgerx484XvNSS/1I4VRbO2KTCHiWax+w9wylhVvvGhZ3O3Cqmc6IuB7eDdi3I7ZigA/O53ducvy3Yq6diEHekVOuvarasccNyyr7y+7firl2ohoGeQetufZM+4Yvxx43LKvsH3NTWms/IebaiWrYatiFtX5+jUOO3i6gs6zS3VB9TSKIrJX5+vNYPSJ7nMl71KnHjRXLKt1lRPDkzk22d0d+Me9O5I5B3iM/eV8vqZ40qywqJo+eRdVlBbr1tWaNO1F3Upeu8XMKlJVdjxu3r/WS6kmzC6Uy8g5pLbM1cOtrDXh//YmoJlV18nYdI7NGJvRFOqdulv20bGkGl6/Fp+LHDORR/HsQJQ3r5Ov8ngIVlG3rRtp20PYzgZPPZSMN8E5pFlbFEIUvVekaL6dAedWa9tm2bgTHz1xsSyXYbYwSAD//jhX45vxbfamn37ZuxPU0JL9uWDKEqwuLnj9/8sGNjmkWVsUQhStVQd6ptNFvhUZr2qdYKjcFUXOj1PRrl/DMy6+3LTAqgHPfL+OJ+zfg8RdO480rFf9/GR8OzZwP5PsYGcHkRzbaBmWnlFQ+l2UgJ4pQqtI1Qe2M9LKjtVyp4ukT844VJBdKZUy/dgmlkAN8bSzeZ91OzAZfTsGau06J4ilVM3m/FTJOvKZ3Oi1pB5lCCVM+l206uNxOUK8tEQUrVUEeCCYHHMSO1vjUNNW4Vd9Yf6m5laAyLUMUP6lK1wRlEHe0Lh82MGzY/3M/vGUUf/ThDY4VP+aahbkWUSyVG03a9kzNYXWHrpxEFB0G+S7Ylf49vGW00fI2bvtbjYxAFbhik5s3hgRjt63Ajs15PLRltG3s1ry63VqEeUfi1pWTiKKTunRNUJxSE3HY+GQlAkCBUtl+gddsMbBjcx5jt63AkW+90aj2yWUNPHbv+sbfs9NahLnngCkbovjgTD5gboGw2xSP21GEnajWArmbC6VyIxVjLedsrYX3UmrKNstE8cIgHzCnQGju5jRTPMuHDeSyRscGZvlcFg9tGbUtT3yoniLqtQnaylzW025gL2sR7ApJFC9M1wTMqR+LWYVil8rw0lNn7LYVjaqW3LABVeDpE/NYmcviyZ2bAKDte3hhjm3P1JztdevM3FomWSyVUc8EtX0vIoqPVDUo65duOl16/Rq3XwjA9QDsxgzOecvzuO1YdaqR77ajJxEFy61BWeqDvNceNHHhJRivGT/iWIefd/g7RdWhk4h65xbkU52u8dKDZs/UHHZPzTkGx35zmqVbH3farOU2K+eOVaJkSl2Qt87ch0RcTycC2uvAAXQd+IJIb2QcxmxdeHVbF3DDHatEyZOqIN86c+8U4FuVK1U8/sLprgK13V1DN780nMZsfZyzciIypSrIe+ke2cmbVyqNWnI/gdqtRNFP8HU7Ms+Ks3IiAlJWJx/GRh2vJ0sFdWAJW/oSkR+pCvJOG3UyIj31oPESqJ2e2+/mIR6ZR0R+pCpd4/fg6NaF0stXF2x7wHgJ1N0uhtphKoaIvEpVkPe7INkaTJ1qyb0Eai6GElEUUr8Zyi/u8iSiuOFmqAAxVUJEgyRVC69ERGnDIE9ElGBM1wSI+XoiipueZvIiMikiZ0TkWyLydyKSs1zbLyKvishZEdne+1Djze6Qa555SkRR6zVd8yKAd6vqewD8B4D9ACAi7wLwMQDrAdwF4E9FpLuz7zwozBaxdeIY1owfwdaJY5EEVi8nKxER9VtPQV5Vv6qqC/UPTwC4tf7n+wB8SVWvqup3AbwK4I5enstJXGbQQbUtICIKUpALr78B4O/rf84DeN1y7Xz9scDFZQYdVNsCIqIgdQzyIvI1Efl3m//us3zOZwEsAHjafMjmW9nuuhKRXSIyLSLTFy9e9P0XiMsMmo3DiCiOOlbXqOoH3K6LyCcA3APg/Xp9++x5AKssn3YrgAsO3/8AgANAbcerhzE3cToFqd8zaLYtIKI46qmEUkTuAvBpAL+oqlcslw4D+BsR+RyAlQDeCeDfenkuJ0E2/uoVd8MSUdz0Wif/BQA3AHhRasfPnVDV31LV0yJyEMC3UUvjfFJVezutwwFn0EREztigjIhowLk1KGNbAyKiBGOQJyJKMAZ5IqIEY5AnIkowBnkiogSLVXWNiFwE8FrU46i7GcD/Rj2IGOHr0Y6vSTO+Hs36+XrcpqojdhdiFeTjRESmnUqS0oivRzu+Js34ejSLy+vBdA0RUYIxyBMRJRiDvLMDUQ8gZvh6tONr0oyvR7NYvB7MyRMRJRhn8kRECcYgT0SUYAzyHojI74mIisjNUY8lSiIyKSJnRORbIvJ3IpKLekxREJG7ROSsiLwqIuNRjydqIrJKRI6LyCsiclpEPhX1mOJARDIiMisiX45yHAzyHYjIKgC/BGA+6rHEwIsA3q2q7wHwHwD2RzyevhORDIAvAvgVAO8C8HEReVe0o4rcAoC9qvozALYA+CRfEwDApwC8EvUgGOQ7exLA78PhjNo0UdWvqupC/cMTqB3rmDZ3AHhVVf9LVa8B+BKA+zp8TaKp6huq+s36n3+IWmBL9ak9InIrgLsB/HnUY2GQdyEi9wIoqurJqMcSQ78B4O+jHkQE8gBet3x8HikPaFYishrAZgAvRzuSyH0etcnhYtQD6fX4v4EnIl8D8JM2lz4L4DMAfrm/I4qW2+uhqs/XP+ezqN2iP93PscWE2DyW+rs8ABCRHwNwCMBuVf1B1OOJiojcA+B7qjojIu+LejypD/Kq+gG7x0VkA4A1AE7Wz6+9FcA3ReQOVf3vPg6xr5xeD5OIfALAPQDer+ncZHEewCrLx7cCuBDRWGJDRAzUAvzTqvpc1OOJ2FYA94rIBwHcCODHReQpVfQrROMAAACxSURBVH04isFwM5RHInIOwJiqprbLnojcBeBzAH5RVS9GPZ4oiMgS1Bad3w+gCOAbAH5VVU9HOrAISW0W9NcALqnq7qjHEyf1mfzvqeo9UY2BOXny4wsA3gbgRRGZE5E/i3pA/VZfeP5tAEdRW2A8mOYAX7cVwK8BuLP+czFXn8VSDHAmT0SUYJzJExElGIM8EVGCMcgTESUYgzwRUYIxyBMRJRiDPBFRgjHIExEl2P8DKO3PcN0auOsAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(x, y)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 损失函数\n",
    "def J(theta, X_b, y):\n",
    "    try: \n",
    "        return np.sum((y - X_b.dot(theta))**2) / len(X_b)\n",
    "    except:\n",
    "        return float(\"inf\")\n",
    "    \n",
    "# 求导   \n",
    "def dJ(theta, X_b, y):\n",
    "    res = np.empty(len(theta))\n",
    "    res[0] = np.sum(X_b.dot(theta) - y)\n",
    "    for i in range(1,len(theta)):\n",
    "        res[i] = (X_b.dot(theta) - y).dot(X_b[:,i]) \n",
    "    \n",
    "    return res * 2 / len(X_b)\n",
    "\n",
    "# 批量梯度下降\n",
    "def gradient_descent(X_b, y, initial_theta, eta, n_iters=1e4, epsilon=1e-8):\n",
    "    theta = initial_theta\n",
    "    i_iter = 0\n",
    "    \n",
    "    while i_iter < n_iters: \n",
    "        gradient = dJ(theta, X_b, y)\n",
    "        last_theta = theta\n",
    "        theta = theta - eta * gradient\n",
    "        \n",
    "        if(abs(J(theta, X_b, y) - J(last_theta, X_b, y)) < epsilon):\n",
    "            break\n",
    "            \n",
    "        i_iter += 1\n",
    "        \n",
    "    return theta\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.92 s, sys: 26.8 ms, total: 1.95 s\n",
      "Wall time: 490 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "X_b = np.hstack([np.ones((len(X),1)), X])\n",
    "initial_theta = np.zeros(X_b.shape[1])\n",
    "eta = 0.01\n",
    "\n",
    "theta = gradient_descent(X_b, y, initial_theta, eta)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.0003629, 4.0022586])"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 随机梯度下降法"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "def dJ_sgd(theta, X_b_i, y_i):\n",
    "    return X_b_i.T.dot(X_b_i.dot(theta) - y_i) * 2. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "def sgd(X_b, y, initial_theta, n_iters):\n",
    "    t0 = 5\n",
    "    t1 = 50\n",
    "    \n",
    "    def learning_rate(t):\n",
    "        return t0 / (t + t1)\n",
    "    \n",
    "    theta = initial_theta\n",
    "    \n",
    "    for cur_iter in range(n_iters):\n",
    "        rand_i = np.random.randint(len(X_b))\n",
    "        gradient = dJ_sgd(theta, X_b[rand_i], y[rand_i])\n",
    "        theta = theta - learning_rate(cur_iter) * gradient\n",
    "        \n",
    "    return theta"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 567 ms, sys: 0 ns, total: 567 ms\n",
      "Wall time: 567 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "X_b = np.hstack([np.ones((len(X),1)), X])\n",
    "initial_theta = np.zeros(X_b.shape[1])\n",
    "\n",
    "theta = sgd(X_b, y, initial_theta, n_iters=len(X_b)//3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([3.00962412, 4.09341491])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "theta"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用我们的SGD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 100000\n",
    "\n",
    "x = np.random.normal(size=m)\n",
    "X = x.reshape(-1,1)\n",
    "y = 4.*x + 3. + np.random.normal(0, 3, size=m)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "from encapsulations.LinearRegression import LinearRegression"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "lin_reg = LinearRegression()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "LinearRegression()"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin_reg.fit_sgd(X, y, n_iters=1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([4.03466837])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin_reg.coef_"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.992087430586617"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "lin_reg.intercept_"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 使用真实数据"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn import datasets\n",
    "boston = datasets.load_boston()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "X = boston.data\n",
    "y = boston.target\n",
    "\n",
    "X = X[y < 50.0]\n",
    "y = y[y < 50.0]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "from encapsulations.model_selection import train_test_split\n",
    "X_train, X_test, y_train, y_test = train_test_split(X, y, seed=666)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.preprocessing import StandardScaler\n",
    "sc = StandardScaler()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "sc.fit(X_train)\n",
    "X_train_standard = sc.transform(X_train)\n",
    "X_test_standard = sc.transform(X_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.71 ms, sys: 4.18 ms, total: 7.89 ms\n",
      "Wall time: 5.84 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.7857275413602652"
      ]
     },
     "execution_count": 23,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from encapsulations.LinearRegression import LinearRegression\n",
    "lin_reg2 = LinearRegression()\n",
    "%time lin_reg2.fit_sgd(X_train_standard, y_train, n_iters=2)\n",
    "lin_reg2.score(X_test_standard,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 111 ms, sys: 4.47 ms, total: 115 ms\n",
      "Wall time: 102 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.808560757055621"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from encapsulations.LinearRegression import LinearRegression\n",
    "lin_reg2 = LinearRegression()\n",
    "%time lin_reg2.fit_sgd(X_train_standard, y_train, n_iters=50)\n",
    "lin_reg2.score(X_test_standard,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 206 ms, sys: 3.99 ms, total: 210 ms\n",
      "Wall time: 203 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.8129434245278827"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from encapsulations.LinearRegression import LinearRegression\n",
    "lin_reg2 = LinearRegression()\n",
    "%time lin_reg2.fit_sgd(X_train_standard, y_train, n_iters=100)\n",
    "lin_reg2.score(X_test_standard,y_test)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## scikit learn 中的SGD"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "from sklearn.linear_model import SGDRegressor\n",
    "sgd_reg = SGDRegressor()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 3.89 ms, sys: 0 ns, total: 3.89 ms\n",
      "Wall time: 2.92 ms\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "SGDRegressor(alpha=0.0001, average=False, early_stopping=False, epsilon=0.1,\n",
       "             eta0=0.01, fit_intercept=True, l1_ratio=0.15,\n",
       "             learning_rate='invscaling', loss='squared_loss', max_iter=1000,\n",
       "             n_iter_no_change=5, penalty='l2', power_t=0.25, random_state=None,\n",
       "             shuffle=True, tol=0.001, validation_fraction=0.1, verbose=0,\n",
       "             warm_start=False)"
      ]
     },
     "execution_count": 27,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "%time sgd_reg.fit(X_train_standard,y_train)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.8122443512755939"
      ]
     },
     "execution_count": 28,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "sgd_reg.score(X_test_standard,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 65 ms, sys: 3.62 ms, total: 68.6 ms\n",
      "Wall time: 67.7 ms\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/birch/anaconda3/lib/python3.7/site-packages/sklearn/linear_model/_stochastic_gradient.py:1187: ConvergenceWarning: Maximum number of iteration reached before convergence. Consider increasing max_iter to improve the fit.\n",
      "  ConvergenceWarning)\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "0.8128784761622174"
      ]
     },
     "execution_count": 29,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from sklearn.linear_model import SGDRegressor\n",
    "sgd_reg2 = SGDRegressor(n_iter_no_change=1e5)\n",
    "%time sgd_reg2.fit(X_train_standard,y_train)\n",
    "sgd_reg2.score(X_test_standard,y_test)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}
